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STRONG  STABILITY  PRESERVING  HIGH-ORDER  TIME  DISCRETIZATION 

METHODS 


SIGAL  GOTTLIEB*,  CHI-WANG  SHUt,  AND  EITAN  TADMOR^ 

Abstract.  In  this  paper  we  review  and  further  develop  a  class  of  strong-stability  preserving  (SSP) 
high-order  time  discretizations  for  semi-discrete  method-of-lines  approximations  of  partial  differential  equa¬ 
tions.  Termed  TVD  (total  variation  diminishing)  time  discretizations  before,  this  class  of  high-order  time 
discretization  methods  preserves  the  strong-stability  properties  of  first-order  Euler  time  stepping  and  has 
proved  very  useful  especially  in  solving  hyperbolic  partial  differential  equations.  The  new  contributions  in 
this  paper  include  the  development  of  optimal  explicit  SSP  linear  Runge-Kutta  methods,  their  application 
to  the  strong  stability  of  coercive  approximations,  a  systematic  study  of  explicit  SSP  multi-step  methods, 
and  a  study  of  the  strong-stability  preserving  property  of  implicit  Runge-Kutta  and  multi-step  methods. 

Key  wrords.  strong-stability  preserving,  Runge-Kutta  methods,  multi-step  methods,  high-order  accu¬ 
racy,  time  discretization 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction,  It  is  a  common  practice  in  solving  time-dependent  Partial  Differential  Equations 
(PDEs)  to  discretize  first  the  spatial  variables  to  obtain  a  semi-discrete  method-of-lines  scheme.  This  is 
then  a  system  of  Ordinary  Differential  Equations  (ODEs)  in  the  time  variable  which  can  be  discretized  and 
solved  by  an  ODE  solver.  A  relevant  question  here  is  stability.  For  problems  with  smooth  solutions,  usually 
a  linear  stability  analysis  is  adequate.  For  problems  with  discontinuous  solutions,  however,  such  as  solutions 
to  hyperbolic  problems,  a  stronger  measure  of  stability  is  usually  required. 

In  this  paper,  we  review  and  further  develop  a  class  of  high-order  strong-stability  preserving  (SSP)  time 
discretization  methods  for  the  semi-discrete  method-of-lines  approximations  of  PDEs.  This  class  of  time 
discretization  methods  was  first  developed  in  [19]  and  [18]  and  was  termed  TVD  (Total  Variation  Diminishing) 
time  discretizations.  It  was  further  developed  in  [6].  The  idea  is  to  assume  that  the  first-order  forward-Euler 
time  discretization  of  the  method-of-lines  ODE  is  strongly  stable  under  a  certain  norm,  when  the  time  step, 
Ai,  is  suitably  restricted,  and  then  try  to  find  a  higher-order  time  discretization  (Runge-Kutta  or  multi- 
step)  that  maintains  strong  stability  for  the  same  norm,  perhaps  under  a  different  time-step  restriction.  In 
[19]  and  [18],  the  relevant  norm  was  the  total  variation  norm:  the  Euler  forward  time  discretization  of  the 
method-of-lines  ODE  was  assumed  TVD,  hence  the  class  of  high-order  time  discretization  developed  there 
was  termed  TVD  time  discretizations.  This  terminology  was  kept  also  in  [6].  In  fact,  the  essence  of  this 
class  of  high-order  time  discretizations  lies  in  its  ability  to  maintain  the  strong  stability  in  the  same  norm 
as  the  first-order  forward  Euler  version,  hence  “strong  stability  preserving” ,  or  SSP,  time  discretization  is  a 
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more  suitable  term  which  will  be  used  in  this  paper. 

We  begin  this  paper  by  discussing  explicit  SSP  methods.  We  first  give,  in  §2,  a  brief  introduction  for 
the  setup  and  basic  properties  of  the  methods.  We  then  move  in  §3  to  our  new  results  on  optimal  SSP 
Runge-Kutta  methods  of  arbitrary  order  of  accuracy  for  linear  ODEs  suitable  for  solving  PDEs  with  linear 
spatial  discretizations.  This  is  used  to  prove  strong  stability  for  a  class  of  well-posed  problems  Ut  =  L(u) 
where  the  operator  L  is  linear  and  coercive,  improving  and  simplifying  the  proofs  for  the  results  in  [13].  We 
review  and  further  develop  the  results  in  [19],  [18]  and  [6]  for  nonlinear  SSP  Runge-Kutta  methods  in  §4  and 
multi-step  methods  in  §5.  Section  6  of  this  paper  contains  our  new  results  on  implicit  SSP  schemes.  It  starts 
with  a  numerical  example  showing  the  necessity  of  preserving  the  strong  stability  property  of  the  method, 
then  it  moves  on  to  the  analysis  of  the  rather  disappointing  negative  results  about  the  non-existence  of  SSP 
implicit  Runge-Kutta  or  multi-step  methods  of  order  higher  than  one.  Concluding  remarks  are  given  in  §7. 

2.  Explicit  SSP  Methods. 

2.1.  Why  SSP  methods?.  Explicit  SSP  methods  were  developed  in  [19]  and  [18]  (termed  TVD  time 
discretizations  there)  to  solve  systems  of  ODEs 

(2.1)  ^u^L(u), 

resulting  from  a  method-of-lines  approximation  of  the  hyperbolic  conservation  law, 

(2.2)  Ut  =  -/(m)i, 

where  the  spatial  derivative,  /(w)a;,  is  discretized  by  a  TVD  finite  difference  or  finite  element  approximation, 
e.g.,  [8],  [16],  [20],  [2],  [9];  consult  [21]  for  a  recent  overview.  Denoted  by  — it  is  assumed  that  the 
spatial  discretization  has  the  property  that  when  it  is  combined  with  the  first-order  forward  Euler  time 
discretization, 

(2.3)  +  AtL{u^), 

then,  for  a  sufficiently  small  time  step  dictated  by  the  CFL  condition, 

(2.4)  At  <  AtpE^ 

the  Total  Variation  (TV)  of  the  one-dimensional  discrete  solution  :=  does  not 

increase  in  time,  i.e.,  the  following,  so  called  TVD  property,  holds 

(2.5)  Ty(u"+i)  <  TV{u^),  TV{u^)  :=  ^  |u;Vi  -  u]\. 

3 

The  objective  of  the  high  order  SSP  Runge-Kutta  or  multi-step  time  discretization  is  to  maintain  the 
strong  stability  property  (2.5)  while  achieving  higher-order  accuracy  in  time,  perhaps  with  a  modified  CFL 
restriction  (measured  here  with  a  CFL  coefficient,  c) 

(2.6)  At<cAtFE- 

In  [6]  we  gave  numerical  evidence  to  show  that  oscillations  may  occur  when  using  a  linearly  stable, 
high-order  method  which  lacks  the  strong  stability  property,  even  if  the  same  spatial  discretization  is  TVD 
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Fig.  2.1.  Second- order  TVD  MUSCL  spatial  discretization.  Solution  after  the  shock  moves  50  mesh  points.  Left:  SSP 
time  discretization;  Right:  non-SSP  time  discretization. 


when  combined  with  the  first-order  forward  Euler  time  discretization.  The  example  is  illustrative,  so  we 
reproduce  it  here.  We  consider  a  scalar  conservation  law,  the  familiar  Burgers’  equation 


with  a  Riemann  initial  data: 


(2.8) 


w(a;,  0) 


1 ,  if  X  <  0 
-0.5,  if  X  >  0. 


The  spatial  discretization  is  obtained  by  a  second-order  MUSCL  [12],  which  is  TVD  for  forward  Euler  time 
discretization  under  suitable  CFL  restriction. 

In  Fig.  2.1,  we  show  the  result  of  using  a  SSP  second-order  Runge-Kutta  method  for  the  time  discretiza¬ 
tion  (left),  and  that  of  using  a  non-SSP  second-order  R.unge-Kutta  method  (right).  We  can  clearly  see  that 
the  non-SSP  result  is  oscillatory  (there  is  an  overshoot). 

This  simple  numerical  example  illustrates  that  it  is  safer  to  use  a  SSP  time  discretization  for  solving 
hyperbolic  problems.  After  all,  they  do  not  increase  the  computational  cost  and  have  the  extra  assurance  of 
provable  stability. 

As  we  have  already  mentioned  above,  the  high-order  SSP  methods  discussed  here  are  not  restricted 
to  preserving  (not  increasing)  the  total  variation.  Our  arguments  below  rely  on  convexity,  hence  these 
properties  hold  for  any  norm.  Consequently,  SSP  methods  have  a  wide  range  of  applicability,  as  they  can 
be  used  to  ensure  stability  in  an  arbitrary  norm,  once  the  forward  Euler  time  discretization  is  shown  to 
be  strongly  stable^  i.e.,  ||u”  -h  AtL{u''^)\\  <  For  linear  examples  we  refer  to  [7],  where  weighted 

SSP  higher-order  discretizations  of  spectral  schemes  are  discussed.  For  nonlinear  scalar  conservation  laws  in 
several  space  dimensions,  the  TVD  property  is  ruled  out  for  high-resolution  schemes;  instead,  strong  stability 
in  the  maximum  norm  is  sought.  Applications  of  L'^-SSP  higher-order  discretization  can  be  found  in  [3], 
[9]  for  discontinuous  Galerkin  and  central  schemes.  Finally,  we  note  that  since  our  arguments  below  are 
based  on  convex  decompositions  of  high-order  methods  in  terms  of  the  first-order  Euler  method,  any  convex 

^By  the  notion  of  strong  stability  we  refer  to  the  fact  that  there  is  no  temporal  growth,  as  opposed  to  the  general  notion  of 
stability  which  allows  a  bounded  temporal  growth,  ||ii”||  <  Const  ■  ||ii^||  with  any  arbitrary  constant,  possibly  Const  >  1. 
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function  will  be  preserved  by  such  high-order  time  discretizations.  In  this  context  we  refer,  for  example,  to 
the  cell  entropy  stability  property  of  high-order  schemes  studied  in  [17],  [15]. 

2.2.  SSP  Runge-Kutta  methods.  In  [19],  a  general  m  stage  Runge-Kutta  method  for  (2.1)  is  written 
in  the  form: 


=  u", 

1-1 

(2.9)  -.Yl  (“».*'“**'*  +  ,  ai,k  >0,  i  = 

k=0 

„n+l 

Clearly,  if  all  the  Pi^s  are  nonnegative,  Pi^k  >  0,  then  since  by  consistency  =  1?  it  follows  that 

the  intermediate  stages  in  (2.9),  amount  to  convex  combinations  of  forward  Euler  operators,  with  At 
replaced  by  ~^At,  We,  thus,  conclude 

Lemma  2.1.  [19].  If  the  forward  Euler  method  (2.3)  is  strongly  stable  under  the  CFL  restriction  (2.4), 
||it” -h  Atl/(u^)||  <  then  the  Runge-Kutta  method  (2.9)  with  Pi^k  >  0  is  SSP,  <  ||w”l|,  provided 

the  following  CFL  restriction  (2.6)  is  fulfilled^ 

(2.10)  At<cAtFE^  c  =  min^hi. 

Pi,k 


If  some  of  the  Pi^kS  are  negative,  we  need  to  introduce  an  associated  operator  L  corresponding  to 
stepping  backward  in  time.  The  requirement  for  L  is  that  it  approximates  the  same  spatial  derivative(s)  as 
L,  but  that  the  strong  stability  property  holds  (“  either  with  respect  to  the  TV  or  another 

relevant  norm),  for  first-order  Euler  scheme,  solved  backward  in  time,  i.e., 

(2.11)  =u^-AtL{u^). 

This  can  be  achieved,  for  hyperbolic  conservation  laws,  by  solving  the  negative  in  time  version  of  (2.2), 

(2.12)  ut  =  f(u)x. 

Numerically,  the  only  difference  is  the  change  of  upwind  direction.  Clearly,  L  can  be  computed  with  the 
same  cost  as  that  of  computing  L.  We  then  have  the  following  lemma. 

Lemma  2.2.  [19].  If  the  forward  Euler  method  combined  with  the  spatial  discretization  L  in  (2.3) 
is  strongly  stable  under  the  CFL  restriction  (2.4),  +  AtL{u^)\\  <  and  if  Euler ^s  method  solved 

backward  in  time  in  combination  with  the  spatial  discretization  L  in  (2.11)  is  also  strongly  stable  under  the 
CFL  restriction  (2.4),  \\u^  —  AtL{u^)\\  <  ||ti”l|,  then  the  Runge-Kutta  method  (2.9)  is  SSP 
under  the  CFL  restriction  (2.6), 

(2.13)  At  <  cAtpE,  c  =  rnin 

\PiM 

provided  pi^kL  is  replaced  by  Pi^kP  whenever  pi^k  is  negative. 

Notice  that,  if  for  the  same  k,  both  and  must  be  computed,  the  cost  as  well  as  storage 

requirement  for  this  k  is  doubled.  For  this  reason,  we  would  like  to  avoid  negative  pi^k  as  much  as  possible. 
However,  as  shown  in  [6],  it  is  not  always  possible  to  avoid  negative  Pi^k^ 
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2.3.  SSP  multi-step  methods.  SSP  multi-step  methods  of  the  form: 

rn 

(2.14)  =  Y.  +  At/?iL(ix"+'-')) ,  ai  >  0, 

were  studied  in  [18].  Since  —  1,  it  follows  that  is  given  by  a  convex  combination  of  forward  Euler 
solvers  with  suitably  scaled  APs,  and  hence,  similar  to  our  discussion  for  Runge-Kutta  methods  we  arrive 
at  the  following  lemma. 

Lemma  2.3.  [18].  If  the  forward  Euler  method  combined  with  the  spatial  discretization  L  in  (2.3) 
is  strongly  stable  under  the  CFL  restriction  (2.4),  ||'u"  +  AtL{u^)\\  <  ||u^||,  and  if  Euler ^s  method  solved 
backward  in  time  in  combination  with  the  spatial  discretization  L  in  (2.11)  is  also  strongly  stable  under  the 
CFL  restriction  (2.4),  -  AtL{u'^)\\  <  then  the  multi-step  method  (2.14)  is  SSP 

under  the  CFL  restriction  (2.6), 

(2.15)  AtKcAtpE,  c  =  min^, 

*  IPm 

provided  PiL{’)  is  replaced  by  AL(-)  whenever  pi  is  negative. 

3.  Linear  SSP  Runge-Kutta  Methods  of  Arbitrary  Order. 

3.1.  SSP  Runge-Kutta  methods  with  optimal  CFL  condition.  In  this  section  we  present  a  class 
of  optimal  (in  the  sense  of  CFL  number)  SSP  Runge-Kutta  methods  of  any  order  for  the  ODE  (2.1)  where 
L  is  linear.  With  a  linear  L  being  realized  as  a  finite  dimensional  matrix  we  denote,  L{u)  =  Lu.  We  will 
first  show  that  the  m-^tage,  m-th  order  SSP  Runge-Kutta  method  can  have,  at  most,  CFL  coefficient  c  =  1 
in  (2.10).  We  then  proceed  to  construct  optimal  SSP  linear  Runge-Kutta  methods. 

Proposition  3.1.  Consider  the  family  of  m- stage,  m-th  order  SSP  Runge-Kutta  methods  (2.9)  with 
nonnegative  coefficients  o>nd  Pi^k  •  I'he  maximum  CFL  coefficient  attainable  for  such  methods  is  the  one 
dictated  by  the  forward  Euler  scheme. 


At  <  AtpE, 


i.e.,  (2.6)  holds  with  maximal  CFL  coefficient  c  =  1. 

Proof.  We  consider  the  special  case  where  L  is  linear,  and  prove  that  even  in  this  special  case  the  maximum 
CFL  coefficient  c  attainable  is  1.  Any  m-stage  method  (2.9),  for  this  linear  case,  can  be  rewritten  as: 


where 


l  +  VA,-,(AiL) 


jfc=0 


.(0) 


i  =  1, ...,  m 


i-l  i-l 

Ai,o  =  Pl,G,  ^  Oti^k^k.Q  +  ^  Pi^k, 

k  =  l  k-0 


i-l  i-l 

^i,k  ~  ^  ^  “  l,...,i  “  1. 

j=k 
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In  particular,  using  induction,  it  is  easy  to  show  that  the  last  two  terms  of  the  final  stage  can  be  expanded 
as 


im,m— 1 


1=1 


k=2 


.l=k+l 


An,m-2  "  ^2  n  f]  Pi, 


k-2 


1=1 


IA~1 


k=l 

For  a  m-stage,  m-th  order  linear  Runge-Kutta  scheme  Am,k  =  (1+1)1  *  Using  Am,m-i  =  Ilzli 
we  can  rewrite 

m  m  /  m  \  /k  —  2 

0^k,k-l 


il=l,l^k 


-  X!  ,  y  ^^Pk,k-2\  JJ  (  11 


l=k+l 


LZ-1 


/=1 


With  the  non-negative  assumption  on  Pi^s  and  the  fact  Am^m-i  =  11^1  A,z-i  ~  hme  Pi^i-i  >  0  for 

all  /.  For  the  CFL  coefficient  c  >  1  we  must  have  >  1  for  all  k.  Clearly,  Am,m-2  =  (^ii),-  is  possible 

under  these  restrictions  only  if  Pk,k-2  =  0  and  =  1  for  all  A;,  in  which  case  the  CFL  coefficient  c  <  1. 


We  remark  that  the  conclusion  of  Proposition  3.1  is  valid  only  if  the  m-stage  Runge-Kutta  method  is 
m-th  order  accurate.  In  [18],  we  constructed  an  m-stage,  first-order  SSP  Runge-Kutta  method  with  a  CFL 
coefficient  c  =  m  which  is  suitable  for  steady  state  calculations. 

The  proof  above  also  suggests  a  construction  for  the  optimal  linear  m-stage,  m-th  order  SSP  Runge- 
Kutta  methods. 

Proposition  3.2.  The  class  ofm  stage  schemes  given  (recursively)  by: 

(3.1)  + 

m-2 

fe=0 

where  ai^o  =  1 

(3.2)  (y.jji^k  ~  )  k  —  l,...,m  2 

^  m—l 

f^m,m—l  ~  j")  ^m,0  =  1  “  ^  ^  0!m,k 

k=l 

is  an  m-order  linear  E.unge-Kutta  method  which  is  SSP  with  CFL  coefficient  c  =  1, 


Proof.  The  first-order  case  is  forward  Euler,  which  is  first-order  accurate,  and  SSP  with  CFL  coefficient 
c  =  1  by  definition.  The  other  schemes  will  be  SSP  with  a  CFL  coefficient  c  =  1  by  construction,  as  long  as 
the  coefficients  are  non-negative. 

We  now  show  that  scheme  (3.1)-(3.2)  is  m-th  order  accurate  when  L  is  linear.  In  this  case  clearly 


/(*) 


(1  +  AiL)Si(°)  =  E 


t.  «(•  -  *)' 


(AtL)' 


,(0) 


i  =  1, m.  —  1, 
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hence  scheme  (3.1)- (3. 2)  results  in 

/  m-2  j  .| 

^(m)  ^  j  ^  ^  +  a;m,m-l 


m! 


j^o  .  A-,=0 


/!;!(m  “  k)\ 


{AtLf  u 


^  7.(0). 


Clearly,  by  (3.2),  the  coefficient  of  {AtL)^  ^  is  coefficient  of  (AtL)^^  is 

m—2 


,  and  the  coefficient  of  {AtLy  is 


1 

+  ^  “m.,-  - 1. 


3=0 


It  remains  to  show  that,  for  1  <  <  m  —  2,  the  coefficient  of  (AtL)^  is  equal  to 


m—2 


(3.3) 


J _ 


j! 


k\(m  —  k)\  '  Au  —  k)\  k\ 

j=k 


This  will  be  shown  by  induction.  Thus,  we  assume  (3.3)  is  true  for  m,  then  for  m  4- 1  we  have,  for  0  <  A)  < 
m-2,  the  coefficient  of  {AtL)^'^^  is  equal  to 

1  .  i!  _  1  /  1  + 

{k  +  l)\{m-k)\  +  “  (fc  +  1)!  [{m-ky.  +  (/ -  fc)!  j 

1  /  1  iV  ^  a  +  1)!) 

(k  +  iy  [{m-ky^  (/  + (/_/,) 

1  /  1  l\  \ 

(fc  +  1)!  + 

1 

-  (fc  +  1)! 

where  in  the  second  equality  we  used  (3.2)  and  in  the  last  equality  we  used  the  induction  hypothesis  (3.3). 
This  finishes  the  proof. 

Finally,  we  show  that  all  the  a’s  are  non-negative.  Clearly  a2,o  =  0:2,1  =  |  >  0.  If  we  assume  amj  >  0 
for  all  j  =  0, ...,  m  —  1,  then 

077i-j-l,j  =  -jCXmJ—l  ^  Oj  j  —  l,...,m  I5  0777-j-1^771 
and,  by  noticing  that  am+ij  <  for  all  j  =  1,  ...,m,  we  have 

m  m 

Om-l-ljO  =  1  —  ^  ^  Om4-l,j  ^  1  ~  ^  ^  —  l  =  0.  B 

j=l  j=l 

As  the  m-stagc,  m-th  order  linear  Runge-Kutta  method  is  unique,  we  have  in  effect  proved  this  unique 
m-stage,  m-th  order  linear  Runge-Kutta  method  is  SSP  under  CFL  coefficient  c  =  1.  If  L  is  nonlinear, 
scheme  (3.1)-(3.2)  is  still  SSP  under  CFL  coefficient  c  =  1,  but  it  is  no  longer  m-th  order  accurate.  Notice 
that  all  but  the  last  stage  of  these  methods  are  simple  forward  Euler  steps. 

We  note  in  passing  the  examples  of  the  ubiquitous  third-  and  forth-order  Runge-Kutta  methods,  which 
admit  the  following  convex  -  and  hence  SSP  decompositions 

3 

(3.4)  Y1  M  ^ 

^=0 

(3.5)  Y.  ^  I  + 
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Table  3.1 


Coefficients  am,j  of  the  SSP  methods  (3.1)-(3,2) 


order  m 

^m,0 

^7n,l 

am, 2 

^m,3 

^m,4 

^771,5 

^m,6  ^m,7 

1 

1 

2 

1 

2 

1 

2 

3 

1 

3 

1 

2 

1 

6 

4 

3 

1 

1 

1 

8 

3 

4 

24 

5 

11 

3 

1 

1 

1 

30 

8 

6 

12 

120 

6 

53 

11 

3 

1 

1 

1 

144 

30 

16 

18 

48 

720 

7  • 

1  103 

53 

11 

3 

1 

1 

1 

280 

144 

60 

48 

72 

240 

5040 

8 

2119 

103 

53 

11 

1 

1 

1  1 

5760 

280 

288 

180 

64 

360 

]440  40320 

We  list,  in  Table  3.1,  the  coefficients  of  these  optimal  methods  in  (3.2)  up  to  m  =  8. 

3.2.  Application  to  coercive  approximations.  We  now  apply  the  optimal  linear  SSP  Runge-Kutta 
methods  to  coercive  approximations.  We  consider  the  linear  system  of  ODEs  of  the  general  form,  with 
possibly  variable,  time-dependent  coefficients, 

(3.6)  ^u{t)  =  L{t)u{t). 

As  an  exami)le  we  refer  to  [7],  where  the  far-from-iiorinal  character  of  the  spectral  differentiation  matrices 
defies  the  straightforward  von-Neumann  stability  analysis  when  augmented  with  high-order  time  discretiza¬ 
tions. 

We  begin  our  stability  study  for  Runge-Kutta  approximations  of  (3.6)  with  the  first-order  forward-Euler 
scheme  (with  {•,  •)  denoting  the  usual  Euclidean  inner  product) 

=  u”  -h  AtnL{t^)u^, 

based  on  variable  time-steps,  :=  Taking  norms  on  both  sides  one  finds 

|^n+l|2  ^  ^2AtnRe{L{nu^,u^)  +  {Atn)^\L{nu^\\ 

and  hence  strong  stability  holds,  \  <  |?z”|,  provided  the  following  restriction  on  the  time  step,  Atn,  is 
met, 

Atn  <  -2Re(L(i^)7/”,7z^)/|L(^^)w^p. 

Following  Levy  and  Tadmor  [13],  we  therefore  make  the 

Assumption  3.1.  (Coercwity) 
r]{t)  >0  such  that 

(3.7) 


The  operator  L{t)  is  (uniformly)  coercive  in  the  sense  that  there  exists 


,  .  .  r-  Re(L(t)u,u) 

T]{t)  :=  mf - >  0. 
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We  conclude  that  for  coercive  L’s,  the  forward  Euler  scheme  is  strongly  stable,  ||/  +  <  1,  if 

and  only  if 


Atn  <  2rj{n. 

In  a  generic  case,  L{t'^)  represents  a  spatial  operator  with  a  coercivity-bound  which  is  proportional 

to  some  power  of  the  smallest  spatial  scale.  In  this  context  the  above  restriction  on  the  time-step  amounts 
to  the  celebrated  Courant-Priedrichs-Levy  (CFL)  stability  condition.  Our  aim  is  to  show  that  the  general 
m-stage,  m-th  order  accurate  Runge-Kutta  scheme  is  strongly  stable  under  the  same  CFL  condition. 


Remark.  Observe  that  the  coercivity  constant,  t],  is  an  upper  bound  in  the  size  of  L;  indeed,  by  Cauchy- 
Schwartz,  ri{t)  <  \L{t)u\  •  |u|/|L(^)?ip  and  hence 

(3.8) 

||L(t)||=supl^|M<^. 

u  |u|  7]{t) 

To  make  one  point 

we  consider  the  fourth-order  Runge-Kutta  approximation  of  (3.6) 

(3.9) 

e  =  L{e)u^ 

(3.10) 

=L(f”+2)(u"  +  ^fc’) 

(3.11) 

z 

(3.12) 

=  L{e+^){u^  +  Atnk^) 

(3.13) 

+  2k^  +  fc-*] . 

Starting  with  second-order  and  higher  the  Runge-Kutta  intermediate  steps  depend  on  the  time  variation 
of  !/(•),  and  hence  we  require  a  minimal  smoothness  in  time,  making 


Assumption  3.2.  (Lipschitz  regidarity) .  We  assume  that  L{-)  is  Lipschitz.  Thus,  there  exists  a  constant 
K  >  0  such  that 


(3.14) 


\\L{t)-Lis)\\<^^\t-s\. 


We  are  now  ready  to  make  our  main  result,  stating 

Proposition  3.3.  Consider  the  coercive  systems  of  ODEs,  (3.6)- (3.7),  with  Lipschitz  continuous  coef¬ 
ficients  (3.14)-  Then  the  fourth-order  B.unge-Kutta  scheme  (3.9-3.13)  is  stable  under  CFL  condition, 

(3.15)  Atn  <  2n{n, 
and  the  following  estimate  holds 

(3.16)  K|  < 


Remark.  The  result  along  these  lines  was  introduced  by  Levy  and  Tadmor  [13,  Main  Theorem],  stating  the 
strong  stability  of  the  constant  coefficients  5-order  Runge-Kutta  scheme  under  CFL  condition 
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Here  we  improve  in  both  simplicity  and  generality.  Thns,  for  example,  the  previous  bound  of  Ci  =  1/31  [13, 
Theorem  3.3]  is  now  improved  to  a  practical  time-step  restriction  with  our  uniform  Cs  =  2. 

Proof.  We  proceed  in  two  steps.  We  first  freeze  the  coefficients  at  t  =  t",  considering  (here  we  abbreviate 

L"  =  L(t")) 


(3.17) 

(3.18) 

(3.19) 

(3.20) 

(3.21) 

Thus, 


/  =  L"u" 


Atn  .1 


A/n  rn\.,n 


f  =  L"(«"  -t-  -^3^)  =  +  -^L^)u 


f  =  U'iu”  + 


J+^in(/+^L") 


=  L”(U"  +  Atnf) 

,,n+l  _ n  . 


P4{AtnL'^)u'^ ,  where  following  (3.5) 


P^iAtnL”)  :=  |/  +  !(/  +  AtL)  -f  i(/  -b  AtL)2  +  ^(/  +  AtL)^ 

Since  the  CFL  condition  (3.15)  implies  the  strong  stability  of  forward-Euler,  i.e.  \\I  -{-  AtnL'^W  <  1,  it  follows 
that  \\P4{AtnL^)\\  <  3/8  +  1/3  +  1/4  +  1/24  =  1.  Thus, 

(3.22)  <  \u^\. 


Next,  we  turn  to  include  the  time  dependence.  We  need  to  measure  the  difference  between  the  exact 
and  the  ‘frozen’  intermediate  values  -  the  k's  and  the  j’s.  We  have 

(3.23)  -j^=0 

(3.24)  e  -  f  =  [L(r+^)  -  L(t’‘)]  (/  +  ^L")u" 

(3.25)  k^-f  =  L(t"+^)^(fc2  -f)  +  [L(t"+i)  -  L(n] 

(3.26)  fc'*  -  f  =  L(t”+^)At„(A;3  -  f)  +  [L(f"+')  -  L{t”)]  Atnf. 


Lipschitz  continuity  (3.14)  and  the  strong  stability  of  forward-Euler  imply 

(3.27) 


Also,  since  ||i"||  <  ,  we  find  from  (3.18)  that  <  |u”|/77(f"),  and  hence  (3.25)  followed  by  (3.27)  and 

the  CFL  condition  (3.15)  imply 


(3.28) 


\k^-f  \  < 


Atn 

27]{t”) 


\k^-f\  + 


K- Atn 

2r/(t") 


Atn 

2ri{V^) 


|u”|  <  <  2K 


Finally,  since  by  (3.19)  f  does  not  exceed,  |/®|  <  ^(pry (1  +  2^(<"))l''^”l’  (3.26)  followed  by  (3.28) 

and  the  CFL  condition  (3.15), 


(3.29) 


At 


7j(f”) 


•/^l  + 


K  ■  Atn 


|7X"|  <  12A'|7/,"|. 
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We  conclude  that  , 


„n+l  ^  ^n+1  ^  ^  ^2{k^  -  f)  +  2{k^  -  f)  +  {k'^  -  f)]  , 

is  upper  bounded  by,  consult  (3.22),  (3.27)-(3.29), 

|yn+l|  <  |i;'‘+l|  +  ^|^2X|M"i+4/i:|u"|  +  12/C|7i"|] 

<  (l  +  3/CAt„)|w"| 


and  the  result  (3.16)  follows.  ■ 

4.  Nonlinear  SSP  Runge-Kutta  Methods.  In  the  previous  section  we  derived  SSP  Runge-Kutta 
methods  for  linear  spatial  discretizations.  As  explained  in  the  introduction,  SSP  methods  are  often  required 
for  nonlinear  spatial  discretizations.  Thus,  most  of  the  research  to  date  has  been  in  the  derivation  of  SSP 
methods  for  nonlinear  spatial  discretizations.  In  [19],  schemes  up  to  third  order  were  found  to  satisfy  the 
conditions  in  Lemma  2.1  with  CFL  coefficient  c  —  1.  In  [6]  it  was  shown  that  all  four  stage,  fourth-order 
Runge-Kutta  methods  with  positive  CFL  coefficient  c  in  (2.13)  must  have  at  least  one  negative  and 
a  method  which  seems  optimal  was  found.  For  large  scale  scientific  computing  in  three  space  dimensions, 
storage  is  usually  a  paramount  consideration.  We  review  the  results  presented  in  [6]  about  strong  stability 
preserving  properties  among  such  low-storage  Runge-Kutta  methods. 

4.1.  Nonlinear  methods  of  second,  third  and  fourth  order.  Here  we  review  the  optimal  (in  the 
sense  of  CFL  coefficient  and  the  cost  incurred  by  L  if  it  appears)  SSP  Runge-Kutta  methods  of  m-stage, 
m-th  order,  for  m  =  2,3,4,  written  in  the  form  (2.9). 

Proposition  4.1,  [6].  If  we  require  >  0,  then  an  optimal  second-order  SSP  Runge-Kutta  method 
(2.9)  is  given  by 

(4.1)  =  u"  +  AtL(u”) 

with  a  CFL  coefficient  c  =  1  in  (2.10).  An  optimal  third-order  SSP  Runge-Kutta  m,ethod  (2.9)  is  given  by 

=  w"  +  AtL(w") 

(4.2)  +  7A<Z/(«^^') 

^  4  4  4 

=  iw"  +  +  ^AtL(M<2)), 

o  O  O 

with  a  CFL  coefficient  c=  1  in  (2.10). 

In  the  fourth-order  case  we  proved  in  [6]  that  we  cannot  avoid  the  appearance  of  negative 

Proposition  4.2.  [6].  The  four-stage^  fourth-order  SSP  Runge  Kutta  scheme  (2.9)  with  a  nonzero 
CFL  coefficient  c  in  (2.13)  must  have  at  least  one  negative 

We  thus  must  settle  for  finding  an  efficient  fourth-order  scheme  containing  L,  which  maximizes  the 
operation  cost  measured  by  where  c  is  the  CFL  coefficient  (2.13)  and  i  is  the  number  of  Ls.  This  way 
we  are  looking  for  a  SSP  method  which  reaches  a  fixed  time  T  with  a  minimal  number  of  evaluations  of  L 
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or  L.  The  best  method  we  could  find  in  [6]  is: 


(4.3) 


with  a  CFL  coefficient  c  =  0.936  in  (2.13).  Notice  that  two  Ls  must  be  computed.  The  effective  CFL 
coefficient,  comparing  with  an  ideal  case  without  Ls,  is  0.936  x  |  =  0.624.  Since  it  is  difficult  to  solve  the 
global  optimization  problem,  we  do  not  claim  that  (4.3)  is  an  optimal  four  stage,  4th-order  SSP  Runge-Kutta 
method. 

4.2.  Low  storage  methods.  For  large  scale  scientific  computing  in  three  space  dimensions,  storage 
is  usually  a  paramount  consideration.  Therefore,  low  storage  Runge-Kutta  methods  [22],  [1],  which  only 
require  two  storage  units  per  ODE  variable,  may  be  desirable.  Here  we  review  the  results  presented  in  [6] 
concerning  strong  stability  preserving  properties  among  such  low-storage  Runge-Kutta  methods. 

The  general  low-storage  Runge-Kutta  schemes  can  be  written  in  the  form  [22],  [1]: 

ri(o)  =  u”,  =  0, 

^  ?:  =  1, . . .  ,m, 

(4.4)  -h  z  =  l,...,m,  By  =  c, 

^n+l 

Only  u  and  du  must  be  stored,  resulting  in  two  storage  units  for  each  variable. 

Following  Carpenter  and  Kennedy  [1],  the  best  SSP  third-order  method  found  by  numerical  search  in 
[6]  is  given  by  the  system 

=  2c^  ^c- 2 

=  36c^  -  36c*^  +  13c2  ~  8c  -f-  4 
=  34c4  -  46c^  +  34c2  -  13c  +  2 

-24(3c-2)(c-l)‘^ 

(3z2  -  zl)^  —  12c(c—  1)(3;22  —  zi) 

+  108(2c  -  l)c''^  -  3(2c  -  1)^5 
242ic(c  -  1)^  -h  72c^g  -b  72c^2c  -  13) 

with  c  =  0.924574,  resulting  in  a  CFL  coefficient  c  -  0.32  in  (2.6).  This  is,  of  course,  less  optimal  than 
(4.2)  in  terms  of  CFL  coefficient,  but  the  low-storage  form  is  useful  for  large  scale  calculations.  Carpenter 
and  Kennedy  [1]  have  also  given  classes  of  five-stage,  fourth-order  low-storage  Runge-Kutta  methods.  We 
have  been  unable  to  find  SSP  methods  in  that  class  with  positive  and  pi^k-  A  low-storage  method  with 
negative  Pi^k  cannot  be  made  SSP,  as  L  cannot  be  used  without  destroying  the  low-storage  property. 


=  Vsec't  +  36^3  _  i35c2  +  84c  - 

12, 

>^2 

^3 

=  12c'‘  -  ISc^*  +  18c2  -  11c  +  2, 

2^4 

=  69c^  -  62c2  +  28c  -  8, 

Z6 

12c(c  -  1)(32:2  -  Zi)  -  {3Z2  - 

21)2 

Bo 

■L>2 

144c(3c-2)(c-l)2 

5 

— 2i  (6c2  —  4c  +  1)  +  323 

Ao 

^2 

(2c+1)2i -3(c  +  2)(2c- 1)2 

=  -h  -  AtL(u^) 


649 


1600 

53989 

2500000 

5121 

~ 20000 


25193600 
102261 


u  - 


5000000 


AtL{u^)  -b 


1600 
4806213 


U)4-§??A^L(u(^)) 


AfL(u/^))  -t- 


23619 

32000' 


(2) 


+ 


20000000 
7873 


7873 
(1) 


10000 


AfL(u/2)) 


+  ^AtLiun  +  ^u(‘)  + 
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4.3.  Hybrid  multi-step  Runge-Kutta  methods.  Hybrid  multi-step  Runge-Kutta  methods  (e.g., 
[10]  and  [14])  arc  methods  which  combine  the  properties  of  Runge-Kutta  and  multi-step  methods.  We 
explore  the  two-step,  two-stage  method: 

(4.5)  —  OL^l'uP'  H-  0:20^”  ^  (/^20^('^^  ^)  +  i^21-^(^”))  j  OL'2h  ^  Oj 

^n+i  _  asoU^~^  -h  031  2  ^  032^” 

(4.6)  +At  (/?3oi(w"“‘)  +  03iL{u^+h  +  h2L{u”))  •  as*  >  0. 

Clearly,  this  method  is  SSP  under  the  CFL  coefficient  (2.10)  if  j3i^k  >  0.  We  could  also  consider  the  case 
allowing  negative  using  instead  (2.13)  for  the  CFL  coefficient  and  replacing  Pi^kL  by  pi^kL  for  the 

negative 

For  third  order  accuracy,  we  have  a  three  parameter  family  (depending  on  c,  030,  and  Q31): 


(4.7) 


(X20  —  3c^  -f-  2c^ 

P20  = 

021  =  1  -  3c2  -  2c^ 


/?21  —  C  +  2c^  -h 

_  2  -f  2030  —  3c  -f-  3o3oc  -|“  03ic^ 

5  -  030  -  30310“^  -  2o3ic^ 

032  =  1  —  031  —  O30 

.  -5  -h  030  4-  9c  -h  3o3oc  -  3o3ic^  -  osic^ 


The  best  method  we  were  able  to  find  is  given  by  c  =  0.4043,  O30  =  0.0605  and  031  =  0.6315,  and  has  a 
CFL  coefficient  c  0.473.  Clearly,  this  is  not  as  good  as  the  optimal  third-order  Runge-Kutta  method  (4.2) 
with  CFL  coefficient  c  =  1.  We  would  hope  that  a  fourth-order  scheme  with  a  large  CFL  coefficient  could 
be  found,  but  unfortunately  this  is  not  the  case  as  is  proven  in  the  following 


Proposition  4.3.  There  are  no  fourth- order  schemes  (4’5)  with  all  non-negative 
Proof.  The  fourth-order  schemes  are  given  by  a  two-parameter  family  depending  on  c,  aao,  and  setting  aai 
in  (4.7)  with 

_  —7  —  0:30  +  10c  —  2a3oc 
”  c2(3  +  8c  +  4c2)  • 

The  requirement  a2i  >  0  enforces,  consult  (4.7),  c  <  -.  The  further  requirement  a2o  >  0  yields 
“I  ^  c  <  |.  0:31  has  a  positive  denominator  and  a  negative  numerator  for  —  ^  <  c  <  and  its  denominator 
is  0  when  c  =  orc  =  -|,  thus,  we  require  -|<c<-|.  In  this  range,  the  denominator  of  asi  is 
negative,  hence  we  also  require  its  numerator  to  be  negative,  which  translates  to  ^  -  Finally,  we 

would  require  032  =  1  —  031  —  0:30  >  0,  which  translates  to  a^o  >  ^  restrictions 

on  030  gives  us  the  following  inequality: 

-7+ 10c  ^  c2(2c  +  l)(2c  +  3)  +  7-10c 
l-h2c  -  (2c-f-l)(2c-l)(c-hl)2  ’ 

which,  in  the  range  of— |<c<— yields  c  >  1  —  a  contradiction.  ■ 
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5.  Linear  and  Nonlinear  Multi-step  Methods.  In  this  section,  we  review  and  further  study  SSP 
explicit  multi-step  methods  (2.14),  which  were  first  developed  in  [18].  These  methods  are  r-th  order  accurate 
if 


(5.1) 


=  1 

i=l 
m 

=  fc  I  j  ,  k  =  l,...,r. 


1=1 


1  =  1 


We  first  prove  a  proposition  which  sets  the  minimum  number  of  steps  in  our  search  for  SSP  multi-step 
methods. 


Proposition  5.1.  For  m>2,  there  is  no  m-step,  order  SSP  method,  and  there  is  no  m-step, 

m-th  order  SSP  method  with  all  non-negative 
Proof.  By  the  accuracy  condition  (5.1),  we  clearly  have 

m  m. 

(5.2)  Ylp{i)ai  =  Y^p'{i)Pi 

1=1  1=1 

for  any  polynomial  p{x)  of  degree  at  most  r  satisfying  p(0)  =  0. 

When  7*  =  m  +  1,  we  could  choose 

(5.3)  p{x)^  j  g{t)dt,  ?(*)  =  !!(*“*)• 

do  i=i 

Clearly  p'(i)  =  q{i)  =  0  for  i  =  1,  We  also  claim  (and  prove  below)  that  all  the  p(i)’s,  i  —  1, . . .  ,m, 

are  positive.  With  this  choice  of  p  in  (5.2),  its  right-hand  side  vanishes,  while  the  left  hand  side  is  strictly 
positive  if  all  >  0  ^ —  a  contradiction. 

When  r  =  m,  we  could  choose 

p{x)  —  x{m  —  x)^~^ . 


Clearly  p{i)  >  0  for  7  =  1,  equality  holds  only  for  i  =  m.  On  the  other  hand,  p'{i)  =  m(l  -  i)(7n  - 
^  equality  holds  only  for  7  =  1  and  i  =  m.  Hence  (5.2)  would  have  a  negative  right  side  and  a 
positive  left  side  and  would  not  be  an  inequality,  if  all  ai  and  pi  are  non-negative,  unless  the  only  nonzero 
entries  are  Pi  and  Pm-  In  this  special  case  we  have  am  =  1  and  =  0  to  get  a  positive  CFL  coefficient 
c  in  (2.15).  The  first  two  order  conditions  in  (5.1)  now  leads  to  Pm  —  m  and  2Pm  —  m,  which  cannot  be 
simultaneously  satisfied. 

We  conclude  with  the  proof  of  the 

Claim.  p{i)  -  q{t)dt  >  0,  q{t)  :=  n”=i  (*  “  i)- 

Indeed,  q{t)  oscillates  between  being  positive  on  the  even  intervals  /o  =  (0,1), /2  =  (2,3),...  and  being 
negative  on  the  odd  intervals,  h  =  (1, 2),  h  =  (3, 4), . . ..  The  positivity  of  the  p(i)’s  for  i  <  (m-f  l)/2  follows 
since  the  integral  of  q(t)  over  each  pair  of  consecutive  intervals  is  positive,  at  least  for  the  first  [(m  4- 1)/2] 
intervals, 

p{2k  +  2)-p{2k)=  [  \q{t)\dt-  [  \q{t)\dt=  f  -  [  \{l  -  t){2  -  t) . . .  {m  -  t)\dt 

—  [  1(1  -  t){2  -  t) . . .  (m  -  1  -  i)|  X  (|(m  -  t)\  -  \t\)dt  >0,  2A:  4- 1  <  (m  -h  l)/2. 

dl2k 
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For  the  remaining  intervals,  we  note  the  symmetry  of  q{t)  w.r.t.  the  midpoint  {m  +  l)/2,  i.e.,  q{t)  = 
+  1  -  t),  which  enables  us  to  write  for  i  >  (m  +  l)/2 


(5.4) 


ri 

p(i)  z=  /  /  q{m-^l-t)dt 

Jo  V(m4-l)/2 

=  /  q{t)dt  +  (-1)'"  /  q{t')dt'. 

Jo  Jm+l—i 


Thus,  if  m  is  odd  then  p(i)  =  p(m.  +  1  “  i)  >  0  for  i>  {m-\-  l)/2.  If  m  is  even,  then  the  second  integral  on 
the  right  of  (5.4)  is  positive  for  odd  Ts,  since  it  starts  with  a  positive  integrand  on  the  even  interval, 

And  finally,  if  m  is  even  and  i  is  odd,  then  the  second  integral  starts  with  a  negative  contribution  from  its 
first  integrand  on  the  odd  interval,  while  the  remaining  terms  that  follow  cancel  in  pairs  as  before;  a 

straightforward  computation  shows  that  this  first  negative  contribution  is  compensated  by  the  positive  gain 
from  the  first  pair,  i.e.,  • 

pQ,  ^77i-f-2 — i 

p{m.  +  2  -  i)  >  /  q{t)dt  4-  /  q{t)dt  >  0,  m  even,  i  odd. 

Jo  Jm+l—i 

This  concludes  the  proof  of  our  claim.  ■ 


We  remark  that  [4]  contains  a  result  which  states  that  there  are  no  linearly  stable  m-step,  (m  4-  l)-th 
order  method  when  m  is  odd.  When  m  is  even,  such  linearly  stable  methods  exist  but  would  require  negative 
tti.  This  is  consistent  with  our  result. 

In  the  remainder  of  this  section  we  will  discuss  optimal  m  step,  m-th  order  SSP  methods  (which  must 
have  negative  pi  according  to  Proposition  5.1  and  m  step,  (m  —  l)“th  order  SSP  methods  with  positive  Pi, 
For  two-step,  second-order  SSP  methods,  a  scheme  was  given  in  [18]  with  a  CFL  coefficient  c  =  | 
(Scheme  1  in  Table  5.1).  We  prove  this  is  optimal  in  terms  of  CFL  coefficients. 

Proposition  5.2.  For  two-step ^  second-order  SSP  methods,  the  optimal  CFL  coefficient  c  in  (2.15)  is 

1 

2  ' 

Proof.  The  accuracy  condition  (5.1)  can  be  explicitly  solved  to  obtain  a  one-parameter  family  of  solutions 

a2  =  1-0^1,  /3i  =  2-“ai,  P2~-^OL\‘ 

The  CFL  coefficient  c  is  a  function  of  and  it  can  be  easily  verified  that  the  maximum  is  c  =  ^  achieved 
at  ai  —  ■ 

We  move  on  to  three-step,  second-order  methods.  It  is  now  possible  to  have  SSP  schemes  with  positive 
ai  and  Pi.  One  such  method  is  given  in  [18]  with  a  CFL  coefficient  c  =  |  (Scheme  2  in  Table  5.1).  We  prove 
this  is  optimal  in  CFL  coefficient  in  the  following  proposition.  We  remark  that  this  multi-step  method  has 
the  same  efficiency  as  the  optimal  two-stage,  second-order  Runge-Kutta  method  (4.1).  This  is  because  there 
is  only  one  L  evaluation  per  time  step  here,  compared  with  two  L  evaluations  in  the  two-stage  Runge-Kutta 
method.  Of  course,  the  storage  requirement  here  is  larger. 


Proposition  5.3.  If  we  require  pi  >  0,  then  the  optimal  three-step,  second-order  method  has  a  CFL 
coefficient  c  =  ^. 
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Proof.  The  coefficients  of  the  three-step,  second-order  method  are  given  by, 

^1  “  2  “  A  +  Pz) )  <^-2  =  -3  -f  2pi  —  2/^3,  03  =  -  (2  -  /?i  -h  /92  4-  3/^3) . 

For  CFL  coefficient  c  >  |  we  need  ^  |  for  all  k.  This  implies 

2q:i  >  ^  6  -  4/^1  -  ^2  +  /^3  >  0 

2q:2  >  ft  =>  —6  -f  4ft  —  ft  —  4ft  >  0 


This  means  that 


ft  “  ft  <  h  —  APi  <  —ft  —  4ft  2ft  <  —3ft. 

Thus,  we  would  have  a  negative  p.  ■ 


We  remark  that  if  more  steps  are  allowed,  then  the  CFL  coefficient  can  be  improved.  Scheme  3  in  Table 

5.1  is  a  four-step,  second-order  method  with  positive  ai  and  ft  and  a  CFL  coefficient  c  =  |. 

We  now  move  to  three-step,  third-order  methods.  In  [18]  we  gave  a  three-step,  third-order  method  with 
a  CFL  coefficient  c  ^  0.274  (Scheme  4  in  Table  5.1).  A  computer  search  gives  a  slightly  better  scheme 
(Scheme  5  in  Table  5.1)  with  a  CFL  coefficient  c  ^  0.287. 

Next  we  move  on  to  four-step,  third-order  methods.  It  is  now  possible  to  have  SSP  schemes  with  positive 
at  and  pi.  One  example  was  given  in  [18]  with  a  CFL  coefficient  c  =  |  (Scheme  6  in  Table  5.1).  We  prove  this 
is  optimal  in  the  CFL  coefficient  in  the  following  proposition.  We  remark  again  that  this  multi-step  method 
has  the  same  efficiency  as  the  optimal  three-stage,  third-order  Runge-Kiitta  method  (4.2).  This  is  because 
there  is  only  one  L  evaluation  per  time  step  here,  compared  with  three  L  evaluations  in  the  three-stage 
Runge-Kutta  method.  Of  course,  the  storage  requirement  here  is  larger. 


Proposition  5.4.  If  we  require  Pi  >  0,  then  the  optimal  four-step,  third-order  method  has  a  CFL 
coefficient  c  =  | . 

Proof.  The  coefficients  of  the  four  step,  third  order  method  are  given  by, 

ai  =  -  (24  —  11/3]  —  2ft  +  ft  —  2ft) ,  q2  =  — 6  4-  3ft  —  -ft  —  ft  +  -ft, 

0:3  =  4  —  -ft  4-  ft  +  -ft  —  3ft,  a4  =  -  (—6  -f  2ft  —  ft  +  2ft  -h  lift) . 

For  a  CFL  coefficient  c  >  |  we  need  ^  |  for  all  k.  This  implies: 

24  -  13ft  -  2ft  4-  ft  -  2ft  >  0,-36  -f  18ft  -  5ft  -  6ft  +  9ft  >  0, 

24  -  9ft  4-  6ft  -{-ft  -  18ft  >  0,-6  -h  2ft  -  ft  -h  2ft  4-  Oft  >  0. 

Combining  these  (9  times  the  first  inequality  plus  8  times  the  second  plus  3  times  the  third)  we  get: 

-40ft  -  36ft  >  0, 


which  implies  a  negative  p.  ■ 


We  again  remark  that  if  more  steps  are  allowed,  the  CFL  coefficient  can  be  improved.  Scheme  7  in  Table 

5.1  is  a  five-step,  third-order  method  with  positive  a,;  and  ft  and  a  CFL  coefficient  c  =  |.  Scheme  8  in  Table 

5.1  is  a  six-step,  third-order  method  with  positive  Oi  and  Pi  and  a  CFL  coefficient  c  =  0.567. 
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Table  5.1 

SSP  multi-step  methods  (2.14) 


# 

steps 

order 

CFL 

A 

m 

r 

c 

1 

2 

2 

1 

2 

4  1 

5’  5 

8  2 

5’  5 

2 

3 

2 

1 

2 

-  0  - 
4)^5  4 

|,o,o 

3 

4 

2 

2 

3 

8  0  0- 

0,0,0 

4 

3 

3 

0.274 

4  2  1 

7’  7  ’  7 

25  20  37 

12  ’  21  ’  84 

5 

3 

3 

0.287 

2973  351  623 

5000’  12.50’  5000 

1297  49  1087 

625  ’  50  ’  2500 

6 

4 

3 

1 

3 

16  0  0  ^ 

16  0  0  - 

7 

5 

3 

1 

2 

^  0  0  0  — 

32  3  ^3  ^3  '-'5  32 

”  0  0  0  “ 

IG  ’  ^3  16 

8 

6 

3  • 

0.567 

108  Q  Q  Q  0  TL 

—  0  0  0  0  — 

9 

4 

4 

0.154 

29  7  1  1 

72  ’  24  ’  4  ’  18 

481  1055  937  197 

192  ’  576  ’  576  ’  576 

10 

4 

4 

0.159 

1989  2893  517  34 

5000’  10000’  2000’  625 

601613  1167  130301  82211 

240000’  640  ’  80000  ’  240000 

11 

6 

4 

0.245 

TIL  0  0  0  ^ 

0  0  0  ^ 

128’  ’  ’  ’128’  8 

12 

5 

4 

0.021 

1557  1  1  2063 

9 

5323561  2659  904987  1567579  n 

32000’  32000’  120’  48000’ 

10 

2304000  ’  2304000  ’  2304000  ’  768000  ’  ^ 

13 

5 

5 

0.077 

117  11 

4  ’  4  ’  24  ’  6  ’  24 

185  851  91  151  199 

64  ’  288  ’24’  96  ’  576 

14 

5 

5 

0.085 

1  13  8  7  3 

4’  50’  25’  50’  100 

52031  26617  1412  14407  6161 

18000’  9000  ’  375  ’  9000  ’  18000 

15 

6 

5 

0.130 

7  3  4  r,  7  1 

20  ’  10  ’  15’  120  ’  40 

291201  198401  88063  n  17969  73061 

108000  ’  86400  ’  43200  ’  43200  ’  432000 

We  now  move  on  to  four-step,  fourth-order  methods.  In  [18]  we  gave  a  four-step,  fourth-order  method 
(Scheme  9  in  Table  5.1)  with  a  CFL  coefficient  c  «  0.154.  A  computer  search  gives  a  slightly  better  scheme 
with  a  CFL  coefficient  c  0.159,  Scheme  10  in  Table  5.1.  If  we  allow  two  more  steps,  we  can  improve  the 
CFL  coefficient  to  c  =  0.245,  Scheme  II  in  Table  5.1. 

Next  we  move  on  to  five-step,  fourth-order  methods.  It  is  now  possible  to  have  SSP  schemes  with  positive 
ai  and  pi.  The  solution  can  be  written  in  the  following  five-parameter  family: 

05  =  1  —  Oi  —  0^2  ~  0:3  —  04,  ~  ^  24/35)  5 


/32  =  ^  (5  -  64ai  -  4502  -  32a3  -  3704  -  96/35) , 
ft  =  ^  (5  +  32qi  +  27o2  +  40o3  +  59q4  +  144ft) , 
ft  =  ^  (55  -  64ai  -  63o2  -  6403  -  55o4  -  96ft) . 
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We  can  clearly  see  that  to  get  ^2  >  0  we  would  need  ai  <  ^,  and  also  /3i  >  ||,  hence  the  CFL  coefficient 
cannot  exceed  c  <  ^  ^  «  0.034.  A  computer  search  gives  a  scheme  (Scheme  12  in  Table  5.1)  with  a 

CFL  coefficient  c  =  0.021.  The  significance  of  this  scheme  is  that  it  disproves  the  belief  that  SSP  schemes 
of  order  four  or  higher  must  have  negative  P  and  hence  must  use  L  (see  Proposition  4.2  for  Runge-Kutta 
methods).  However,  the  CFL  coefficient  here  is  probably  too  small  for  the  scheme  to  be  of  much  practical 
use. 

We  finally  look  at  five-step,  fifth-order  methods.  In  [18]  a  scheme  with  CFL  coefficient  c  —  0.077  is 
given  (Scheme  13  in  Table  5.1).  A  computer  search  gives  us  a  scheme  with  a  slightly  better  CFL  coefficient 
c  ^  0.085,  Scheme  14  in  Table  5.1.  Finally,  by  increasing  one  more  step,  one  could  get  [18]  a  scheme  with 
CFL  coefficient  c  =  0.130,  Scheme  15  in  Table  5.1. 

We  list  in  Table  5.1  the  multi-step  methods  studied  in  this  section. 

6.  Implicit  SSP  Methods. 

6,1.  Implicit  TVD  stable  scheme.  Implicit  methods  are  useful  in  that  they  typically  eliminate  the 
step-size  restriction  (CFL)  associated  with  stability  analysis.  For  many  applications,  the  backward-Euler 
method  possesses  strong  stability  properties  that  we  would  like  to  preserve  in  higher-order  methods.  For 
example,  it  is  easy  to  show  a  version  of  Harten’s  lemma  [8]  for  the  TVD  property  of  the  implicit  backward- 
Euler  method: 

Lemma  6.1.  (Harten).  The  following  implicit  backward-Euler  method 

(6.1) 

where  Cjj^i  and  are  functions  ofu^  and/or  at  various  (usually  neighboring)  grid  points  satisfying 

(6.2)  >  0,  >  0, 

is  TVD  in  the  sense  of  (2.5)  for  arbitrary  At. 

Proof.  Taking  a  spatial  forward  difference  in  (6.1)  and  moving  terms,  one  gets 

[1  +  At  -  nj+i) 

=  -  ^Iti)  +  {n’r  -  ■ 

Using  the  positivity  of  C  and  D  in  (6.2),  one  gets 

which,  upon  summing  over  would  yield  the  TVD  property  (2.5).  ■ 

Another  example  is  the  cell  entropy  inequality  for  the  square  entropy,  satisfied  by  the  discontinuous 
Galerkin  method  of  arbitrary  order  of  accuracy  in  any  space  dimensions  when  the  time  discretization  is  by 
a  class  of  implicit  time  discretization  including  backward-Euler  and  Crank-Nicholson,  again  without  any 
restriction  on  the  time  step  At  [11]. 

As  in  Section  2  for  explicit  methods,  here  we  would  like  to  discuss  the  possibility  of  designing  higher- 
order  implicit  methods  that  share  the  strong  stability  properties  of  backward-Euler  without  any  restriction 
on  the  time  step  At. 
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Unfortunately,  we  are  not  as  lucky  in  the  implicit  case.  Let  us  look  at  a  simple  example  of  second-order 
implicit  Runge-Kutta  methods: 


(6.3)  yU)  ~  u'^ 

=  a2,ow^  +  -h  /?2 ). 


Notice  that  we  have  only  a  single  implicit  L  term  for  each  stage  and  no  explicit  L  terms,  in  order  to  avoid 
time  step  restrictions  necessitated  by  the  strong  stability  of  explicit  schemes.  However,  since  the  explicit 
term  is  contained  indirectly  in  the  second  stage  through  the  term,  we  do  not  lose  generality  in 
wiitiiig  the  schemes  as  the  form  in  (6.3)  except  for  the  absence  of  the  terms  in  both  stages. 

To  simplify  our  example  we  assume  L  is  linear.  Second-order  accuracy  requires  the  coefficients  in  (6.3) 
to  satisfy 


(6.4) 


<^2,1 


1 

2/3i(l-A)’ 


0^2,0  —  1  —  0:2,1, 


02  = 


1-2A 

2(1  -  A)  ‘ 


To  obtain  a  SSP  scheme  out  of  (6.4)  we  would  require  0:2,0  and  a2,i  to  be  non-negative.  We  can  clearly  see 
that  this  is  impossible  as  0:2,1  is  in  the  range  [4, -foo)  or  (— oo,0). 

We  will  use  the  following  simple  numerical  example  to  demonstrate  that  a  non-SSP  implicit  method 
may  destroy  the  non-oscillatory  property  of  the  back  ward- Euler  method,  despite  the  same  underlying  non- 
oscillatory  spatial  discretization.  We  solve  the  simple  linear  wave  equation 


(6.5) 


Ut  -  Ux 


with  a  step-function  initial  condition: 


(6.6) 


u{x,Q)  — 


1,  if  X  <  0 
0,  if  X  >  0. 


Ux  in  (6.5)  is  approximated  by  the  simple  first  order  upwind  difference: 


(Wj+l  -  Uj) 


The  backward-Euler  time  discretization 


,„n+l  ^  ^  AtL(M"+l) 

for  this  problem  is  unconditionally  TVD  according  to  Lemma  6.1.  We  can  see  on  the  left  of  Fig.  6.1  that  the 
solution  is  monotone.  However,  if  we  use  (6.3)-(6.4)  with  A  =  2  (which  results  in  positive  02  =  |,  02,0  =  f , 
but  a  negative  a2,i  =  -  j)  as  the  time  discretization,  we  can  see  on  the  right  of  Fig.  6.1  that  the  solution  is 
oscillatory. 

In  the  next  two  subsections  we  discuss  the  rather  disappointing  negative  results  about  the  non-existence 
of  high  order  SSP  Runge-Kutta  or  multi-step  methods. 

6.2.  Implicit  Runge-Kutta  methods.  A  general  implicit  Runge-Kutta  method  for  (2.1)  can  be 
written  in  the  form 

=  7/", 

i~l 

(6-7)  7r«  =  q;,*  >0,  i  = 

k=0 

yn+1 
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Fig.  6.1.  First-order  upwind  spatial  discretization.  Solution  after  100  time  steps  at  CFL  number  ^  =  1.4.  Left:  first- 
order  backward- Euler  time  discretization;  Right:  non-SSP  second-order  implicit  Runge-Kutta  time  discretization  (6.3)-(6.4) 
with  fii  =  2. 


Notice  that  we  have  only  a  single  implicit  L  term  for  each  stage  and  no  explicit  L  terms.  This  is  to  avoid 
time  step  restrictions  for  strong  stability  properties  of  explicit  schemes.  However,  since  explicit  L  terms  are 
contained  indirectly  beginning  at  the  second  stage  from  u  of  the  previous  stages,  we  do  not  lose  generality 
in  writing  the  schemes  as  the  form  in  (6.7)  except  for  the  absence  of  the  terms  in  all  stages.  If  these 

L{u^^^)  terms  are  included,  we  would  be  able  to  obtain  SSP  Runge-Kutta  methods  under  restrictions  on 
similar  to  explicit  methods. 

Clearly,  if  we  assume  that  the  first-order  implicit  Euler  discretization 

(6.8)  + A^L(u”+^) 

is  unconditionally  strongly  stable,  ||n”+^||  <  ||w'^^|i,  then  (6.7)  would  be  unconditionally  strongly  stable  under 
the  same  norm  provided  A  >  0  for  all  i.  If  Pi  becomes  negative,  (6.7)  would  still  be  unconditionally  strongly 
stable  under  the  same  norm  if  piL  is  replaced  by  PiL  whenever  the  coefficient  Pi  <  0,  with  L  approximates  the 
same  spatial  derivative(s)  as  L,  but  is  unconditionally  strongly  stable  for  first-order  implicit  Euler,  backward 
in  time: 

(6.9)  AtZ(n^+'). 

As  before,  this  can  again  be  achieved,  for  hyperbolic  conservation  laws,  by  solving  (2.12),  the  negative  in 
time  version  of  (2.2).  Numerically,  the  only  difference  is  the  change  of  upwind  direction. 

Unfortunately,  we  have  the  following  negative  result  which  completely  rules  out  the  existence  of  SSP 
implicit  Runge-Kutta  schemes  (6.7)  of  order  higher  than  one. 


Proposition  6.1.  If  (6.7)  is  at  least  second-order  accurate^  then  ai^k  cannot  be  all  non-negative. 
Proof.  We  prove  that  the  statement  holds  even  if  L  is  linear.  In  this  case  second-order  accuracy  implies 


(6.10) 


i-1 

'y  ^  “  1)  “  1 5 

k~0 


Y 

m 


1 

2 
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where  Xm  and  Ym  can  be  recursively  defined  as 


m— 1 


m-1 


(6.11)  —  /?i,  1^1  —  Pii  Xjn  —  Pm  d"  ^  ^  ^7i  —  Pm^m  d"  ^  ^  ^m,iY^i' 


i=l 


i=l 


Wc  now  show  that,  if  ai^k  ^  0  for  all  i  and  then 


Xm  Ym 


(6.12)  ^^m  -^m  ^  2? 

which  is  clearly  a  contradiction  to  (6.10).  In  fact,  we  use  induction  on  m  to  prove 

(6.13)  (1  -  a)Xm  -  Ym  <  Cm{l  “  for  any  real  number  a, 


where 

(6.14)  c,  =  Ci+i  = 

It  is  easy  to  show  that  (6.14)  implies 

(6-15)  —  =  Cl  <  C2  <!  •  •  •  <  Cm  2‘ 

We  start  with  the  case  m  =  1.  Clearly, 

(1  -  a)Xi  -  Fi  =  (1  -  a)A  -Pi<  ^(1  -  af  =  c,(l  -  af 
for  any  a.  Now  assume  (6.13)-(6.14),  lienee  also  (6.15),  is  valid  for  all  m  <  k,  for  m  =  k  we  have 


k-\ 


(1  -  a)Xk  -  Fc  =  (1  -  a  -  Pk)l3k  +  ^  ock,i  [(1  -  a  -  (ik)Xi  -  F] 

i-l 

<  (1  -  a  -  fik)Pk  +  Ck-i{l  -a-  fikf' 


< 


4(1-C*:_l) 

=  Cfc(l  -a)^ 


[l-af 


where  in  the  first  equality  we  used  (6.11),  in  the  second  inequality  we  used  (6.10)  and  the  induction  hypothesis 
(6.13)  and  (6.15),  and  the  third  inequality  is  a  simple  maximum  of  a  quadratic  function  in  pk^  This  finishes 
the  proof.  ■ 


We  remark  that  the  proof  of  Proposition  6.1  can  be  simplified,  using  existing  ODE  results  in  [5],  if  all 
PiS  are  non-negative  or  all  PiS  are  non-positive.  However,  the  case  containing  both  positive  and  negative 
PiS  cannot  be  handled  by  existing  ODE  results,  as  L  and  L  do  not  belong  to  the  same  ODE. 

6.3.  Implicit  multi-step  methods.  For  our  purpose,  a  general  implicit  multi-step  method  for  (2.1) 
can  be  written  in  the  form 

m 

(6.16)  =  X)  +  At^oI'(«"+^),  Qi  >  0. 

i=l 

Notice  that  we  have  only  a  single  implicit  L  term  and  no  explicit  L  terms.  This  is  to  avoid  time  step 
restrictions  for  norm  properties  of  explicit  schemes.  If  explicit  L  terms  are  included,  we  would  be  able  to 
obtain  SSP  multi-step  methods  under  restrictions  on  At  similar  to  explicit  methods. 


2i 


Clearly,  if  we  assume  that  the  first-order  implicit  Euler  discretization  (6.8)  is  unconditionally  strongly 
stable  under  a  certain  norm,  then  (6.16)  would  be  unconditionally  strongly  stable  under  the  same  norm 
provided  that  /3o  >  0.  If  is  negative,  (6.16)  would  still  be  unconditionally  strongly  stable  under  the  same 
norm  if  L  is  replaced  by  L. 

Unfortunately,  we  have  the  following  negative  result  which  completely  rules  out  the  existence  of  SSP 
implicit  multi-step  schemes  (6.16)  of  order  higher  than  one. 

Proposition  6.2.  If  (6,16)  is  at  least  second-order  accurate,  then  ai  cannot  he  all  non-negative. 
Proof.  Second-order  accuracy  implies 

m  m  m 

(6.17)  '^ai  =  l,  Y^iai  =  0u,  =  0. 

2=1  t=l  2=1 

The  last  equality  in  (6.17)  implies  that  Oj  cannot  be  all  non-negative.  ■ 

7.  Concluding  Remarks.  We  have  systematically  studied  strong  stability  preserving,  or  SSP,  time 
discretization  methods,  which  preserve  strong  stability  of  the  forward- Euler  (for  explicit  methods)  or  the 
backward-Euler  (for  implicit  methods)  first-order  time  discretizations.  Runge-Kutta  and  multi-step  methods 
are  both  investigated.  The  methods  listed  here  can  be  used  for  method-of-lines  numerical  schemes  for  partial 
differential  equations,  especially  for  hyperbolic  problems. 
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